#include <QtCore/QCoreApplication>
//#include <matrix.c>
//extern "C" {
    #include <kfc.h>
//}
#include <iostream>

#define RUN_LU_TESTS false

#define MAT_SIZE 4

void testLUSolve4xN()
{
    int rows = 4, cols = 3;
    srand((int)time(NULL));

    std::cout << "\nMatrix 'A':\n";
    MTX A = alloc_MTX(rows,rows);

    A[0][0] = 45.28712988; A[0][1] = 45.49704797; A[0][2] = 30.05309633; A[0][3] = 45.65405874;
    A[1][0] = 21.54842841; A[1][1] = 97.79320364; A[1][2] = 3.834773928; A[1][3] = 62.46485762;
    A[2][0] = 98.98824078; A[2][1] = 70.35034772; A[2][2] = 42.19445888; A[2][3] = 92.41798687;
    A[3][0] = 54.22910615; A[3][1] = 23.31723987; A[3][2] = 80.16678345; A[3][3] = 35.32604231;
    printmat(rows, rows, A);

    std::cout << "\nMatrix 'B':\n";
    MTX B = alloc_MTX(rows,cols);

    B[0][0] = 70.8916239;  B[0][1] = 28.801952;   B[0][2] = 70.6034163;
    B[1][0] = 50.28415814; B[1][1] = 9.812234578; B[1][2] = 79.79188494;
    B[2][0] = 98.29242652; B[2][1] = 62.0138899;  B[2][2] = 3.164443113;
    B[3][0] = 82.74165666; B[3][1] = 9.859040002; B[3][2] = 96.89635716;
    printmat(rows, cols, B);

    std::cout << "\nResult:\n";
    MTX res = alloc_MTX(rows,cols);
    lu_solve_4_n(A, cols, B, res);
    printmat(rows, cols, res);

    std::cout << "\nDesired result:\n";
    MTX C = alloc_MTX(rows,cols);

    C[0][0] = 5.141590563;  C[0][1] = 2.904538349;  C[0][2] = 9.553361041;
    C[1][0] = 4.133366044;  C[1][1] = 1.554252649;  C[1][2] = 9.97598788;
    C[2][0] = -0.380051347; C[2][1] = -0.872932507; C[2][2] = -0.394420805;
    C[3][0] = -7.416437848; C[3][1] = -3.224594705; C[3][2] = -17.61213221;
    printmat(rows, cols, C);

    dealloc_MTX(rows, A);
    dealloc_MTX(rows, B);
    dealloc_MTX(rows, res);
    dealloc_MTX(rows, C);
}

void testMatInverse()
{
    srand((int)time(NULL));

    std::cout << "\nMatrix 'A':\n";
    MTX A = alloc_MTX(MAT_SIZE,MAT_SIZE);
    matfill(MAT_SIZE,MAT_SIZE,A);
    printmat(MAT_SIZE, MAT_SIZE, A);

    std::cout << "\nVector 'z':\n";
    VECT z = alloc_VECT(MAT_SIZE);
    vectfill(MAT_SIZE, z);
    printvect(MAT_SIZE, z);

    std::cout << "\nVector 'y = A^(-1) * z':\n";
    VECT y = alloc_VECT(MAT_SIZE);
    lu_solve_4_1(A, z, y);
    printvect(MAT_SIZE, y);

    dealloc_MTX(MAT_SIZE, A);
    free(z); free(y);
}

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

#if RUN_LU_TESTS
    std::cout << "\n\n=== Testing 'LU_SOLV_4_1' ===\n";
    testMatInverse();

    std::cout << "\n\n=== Testing 'LU_SOLV_4_N' ===\n";
    testLUSolve4xN();
#endif

    a.exit(0);
}
